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AN  ANALYSTS  OF  AN  IMPLICIT  FACTORED  SCHEME 
FOR  SIMULATING  SHOCK  WAVES 


M.B.  Tyndall 


ABSTRACT 


A  detailed  derivation  and  analysis  of  an  implicit  factored  scheme  (Eeom  and 
Warming  1976)  is  given.  A  one  dimensional  shock  tube  problem  is  solved  numerically 
using  the  factored  scheme.  The  results  and  exact  solution  are  presented  for  this 
problem.  An  analysis  of  the  features  of  the  method  i3  made  and  the  limitations  of  this 
implicit  factored  scheme  for  more  general  applications  to  shock  waves  in  solids  are 
discussed.  An  alternative  approach  which  has  been  developed  by  Boris  and  Book  (1976) 
appears  to  have  wider  applicability  than  the  method  studied  here. 
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the  Jacobian  matrix  of  the  flux  vector  F 
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the  identity  matrix 

subscript  denoting  the  grid  points  in  the  x-direction 
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matrix  of  eigenvectors  of  A 

time 
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the  shock  velocity 
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the  ratio  of  specific  heat* 

a  backward  difference  operator 

a  forward  difference  operator 

the  time  increment 

the  space  increment 

the  classical  forward  deference  operator 

total  energy  per  unit  mass 

parameter  related  to  the  choice  cf  time  differencing 
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parameter  related  to  the  choice  of  time  differencing 
the  density 

the  density  in  region  i 

the  classical  backward  difference  ope.ator 
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1,  INTRODUCTION 


In  this  report  a  detailed  analysis  of  an  implicit  finite-difference  algorithm  (Beam  and 
Warming  11*76)  for  solving  nonlinear  hyperbolic  systems  in  conservative  form  is  giver. 
This  method  is  one  of  many  finite-difference  schemes  used  to  solve  systems  of  hyperbolic 
equations.  Familiar  schemes  for  solving  these  systems  include  the  finite-difference  methods 
of  Godunov,  Lax  and  Wendroff,  McCormack,  the  upwind  scheme,  the  hybrid  scheme  of 
Harten  and  Zwaa,  and  the  flux  corrected  transport  (FCT)  method  of  Boris  and  Book.  An 
excellent  review  of  these  methods  and  others  is  given  by  Sod  1978.  The  cubic-iruerpolated 
pseudo  particle  (CIP)  method  (Takewaki  and  Yabe  1987)  and  finite  element  flux  corrected 
transport  (Limner  and  others  1987)  are  just  two  of  the  more  recent  methods  in  use.  All 
of  these  methods  arc  moat  commonly  used  to  solve  the  Euler  and  Navier-Stokes  equations. 
This  report  looks  at  a  one-dimensional  shock  tube  problem  using  the  Eulerian  (in viscid) 
gasdynamic  equations. 

This  ana'ysis  was  prompted  by  the  need  to  find  a  suitable  numerical  method  that 
could  efficiently  and  accurately  model  shock  waves  in  two-dimensional  solids  Beam  and 
Warming  1976  presented  and  solved  problems  involving  shock  waves  in  two-dimensional 
inviscid  gases,  such  ns  transonic  aerodynamics  and  shock  boundary-layer  interactions.  For 
this  report  a  onc  uir  »ensi'»nat  form  of  tim  algorithm  given  in  Seam  and  Warming  1970  was 
used  to  model  shock  waves  in  an  inviscid  gas.  Based  on  the  performance  in  this  problem 
this  numerical  method  was  assessed  for  its  suitability  to  model  the  more  difficult  problem  of 
shocks  in  two-dimensional  solids,  where  the  solids  arc  assumed  to  behave  hydrodyn&micaily. 

There  are  numerous  difficulties  that  are  encountered  when  trying  to  model  shock  waves 
numerically.  Some  are  related  to  the  specific  problem,  eg.  sleep  gradients  and  boundary 
conditions,  others  to  the  numerical  method  that  is  being  used,  eg.  damping,  dispersion  and 
non-pliycical  oscillation*  (Gibbs  err^r)  while  others  effect  the  whole  system  eg.  accuracy 
and  stability  The  main  aim  of  choosing  a  particular  numerical  method  is  l h it t  it  must 
minimize  these  difficulties.  In  choosing  a  suitable  finite-difference  algorithm  there  is  a  set 
of  bas.c  requirements  to  be  satisfied.  A  complete  list  is  given  by  Boris  and  Book  1976b 
The  requirements  of  most  importance  are: 

1.  Uxi»ct  c.' •loervrttjon  properties  of  the  physical  ••{(nation*  should  be  mirrored  in  the  nu¬ 
merical  method  used. 

2  1  he  numerical  method  must  be  stable  for  a  range  of  v;rid  spacing*  and  lime  steps. 

T  The  density  p  should  remain  positive. 
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4.  The  numerical  method  should  provide  at  least  second  order  accuracy  in  regions  of  the 
problem,  eg.  the  shock  front. 

5.  The  numerical  method  should  be  fast  and  efficient. 

These  requirements  are  listed  in  order  of  decreasing  importance.  The  first  three  require¬ 
ments  must  be  satisfied  by  the  fir.ite-difference  algorithm  otherwise  the  results  cannot  be 
considered  to  be  realistic.  As  well  as  these  requirements  the  more  difficult  problems  will 
require  a  finite-difference  algorithm  that  is  fairly  robust  and  can  cope  with  a  wide  range 
of  boundary  conditions  and  equations  of  state.  For  modelling  shock  waves  in  solids, 
elastic-plastic  terms  must  also  be  included  and  tolerated  by  the  algorithm. 

An  implicit  finite  difference  scheme  has  been  chosen  here  to  solve  the  one-dimensional 
form  of  the  inviscid  gasdynamic  equations.  Long  computational  times  have  excluded  the 
possibility  of  using  an  explicit  finite  difference  scheme  since  the  stability  bound  of  an 
explicit  algorithm  forces  a  time  step  that  can  be  orders  of  magnitude  smaller  ihan  that 
required  for  accuracy.  There  is  a  limited  number  of  spatial  difference  approximations 
that  can  be  used  for  the  conservative  form  of  the  inviscid  gasdynamic  equations.  Only 
centered  difference  operators  lead  to  difference  methods  that  are  simultaneously  stable  for 
both  positive  and  negative  characteristic  speeds  (i.e.  eigenvalues)  that  are  associated  with 
spatial  flux  terms.  One-sided  (or  upwind)  schemes  can  be  used  when  the  flux  Verms  are  split 
into  components  corresponding  to  either  their  negative  or  positive  characteristic  speeds. 
One-sided  schemes  have  superior  dissipative  and  dispersive  properties  compared  to  those 
of  centered  schemes  (Steger  and  Warming  1981).  Therefore  a  more  robust  and  efficient 
algorithm  can  be  obtained  by  splitting  the  flux  terms  and  applying  one-sided  differences 
The  implicit  factored  scheme  developed  by  Ream  and  Warming  1976  incorporater  these 
principles  and  it  the  more  commonly  used  method  of  this  type. 

In  the  second  section  a  brief  description  nf  the  implicit  factored  scheme  of  Beam  and 
Warming  1976  is  given.  The  Iticmann  problem  in  one  dimension  is  solved  numerically  with 
this  method.  A  discussion  of  this  problem  is  given  in  Section  3  and  the  results  are  presented 
in  the  Section  4.  Section  5  sums  up  the  analysis  end  compares  'his  method  with  others 
and  conclusions  are  drawn  concerning  the  suitability  of  the  Beam  and  Warming  algorithm 
for  modelling  shock  waves  in  solids. 
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In  this  section  a  brie?  summary  ol  Ihe  one- dimension  a}  form  of  the  method  in  Deem 
and  Warming  1976  is  given.  A  full  derivation  is  presented  in  Appendix  A. 

2.1  One  Dimension 

In  one  dimension  the  Eulrrian  invisrid  gasdyr.amic  equations  (Appendix  A),  namely 
the  continuity  equation,  conservation  of  mass  and  conservation  of  energy,  c  >  be  writter 
in  conservation  (or  vector)  form  as 


dU  9F{V) 
8t  8z 


(1) 


(  p\  (  Eu  \ 

U  =  I  /m  F(U)  =  pu1  +  p  I 

\e)  V(£  +  P)u/ 


(2) 


p  is  the  density,  u  the  particle  velocity,  E  the  total  energy  per  unit  volume  and  the  pressure 
p  =  (7  -  l)l£  -  Jpu2),  7  is  the  ratio  of  specific  heats.  An  ideal  gas  equation  of  state  is 
assumed. 

Equation  (1)  is  the  futiu„mental  system  to  solve  It  is  a  form  of  the  Enin  equation  for 
non  viscous  fluids  The  scheme  rlivised  by  Beam  and  Warming  1976  uses  a  factorisation  of 
the  equation  after  time  differencing.  This  approach  relies  on  the  explicit  form  of  equation 
(1)  including  the  use  of  an  equation  of  state  that  can  he  expressed  in  the  functional  form 

P  =  p/(<)  <3) 

where  e  is  the  internal  energy  per  unit  mass  and  /  is  a  function.  When  the  equation  of 
s'ste  can  he  expressed  in  this  functional  form  then  the  nonlinear  flux  vector  F(U)  is  a 
homogeneous  func'ion  of  degree  one  Using  this  homogeneity  property 

F(U)=/t(U)U  (4) 


of  the  Euler  equations,  where  A  is  the  Jacobian  matrix  of  F,  a  linear  tir.ie-difiVreii'-ed  form 
of  equation  (1)  can  be  w-itten  (Appendix  A)  as 
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where  I  is  the  identity  matrix,  AUn  =  Un  +  I  -  Un,  At  the  discrete  time  increment  and 
U(t)  =  U(nAt)  =  Un.  The  parameter*  8  and  (  define  the  particular  time-difference 
approximation  that  can  be  used.  Example*  of  some  familiar  implicit  schemes  are 
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trapeioidal  formula; 

9  =  1, 

o 
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backward  Euler; 

9  =  1, 

three-point  backward 

The  Jacobian  matrix  A  has  the  eigenvalues  u  and  vie  where  c  =  ( "T )  is  the  lo¬ 
cal  6pc*ed  of  sound.  One-sided  spatial  differences  have  superior  properties  when  compared 
lo  centered  rpaltal  difference >  (Section  1).  For  this  reason  Beam  and  Warming  1976  use 
one-sided  difference*  in  their  scheme.  To  be  able  to  apply  these  o.ie-sided  spatial  differ¬ 
ences  the  flux  terms  must  be  split  according  to  the  sign  of  their  characteristic  speeds  (i.e. 
eigenvalues).  The  flux-vector  F  can  be  split  into  two  parts,  F+  and  F-,  (6ee  Appendix  A) 
F*  corresponds  to  the  positive  eigenvalues  of  A  and  F“  to  the  negative  eigenvalues  of  A 
The  partial  derivative  •!-  in  equation  (5)  can  be  replaced  with  one-sided  spatial  difference 
appro cimationa.  To  maintain  the  stability  of  the  system  a  backward  difference  operator  is 
used  for  positive  eigenvalues  of  A  and  a  forward  difference  operator  for  negative  eigenvalues 
of  A.  By  dropping  the  third  order  term  0(At*),  splitting  the  flux  vectors  and  applying 
one-sided  spatial  differences  to  equation  (5),  a  noniterative  implicit  finite  difference  «-chcme 
for  equation  (1),  that  is  second  order  accurate,  can  be  written  as. 
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where  *  -  j'Av,  V*  is  the  d  j'sical  backward  difference  operator  defined  in  Appendix  A 
(equation  (/t2ti))  and  Ax  the  Us'ioil  forward  difference  operator  defined  in  Appendix  A 
(equation  (/I27j);  .fj  and  &*  defined  in  Appendix  A  (after  conation  (.128))  The  flow 
vector  U  in  equation  (1)  has  been  replaced  bv  U*‘,  where  the  superscript  n  d«-n«-t»s  the 
time  level  and  the  subscript  j  denotes  the  grid  point  location.  The  'joiatnm  -i  ;n  means 
the  value  of  A  evaluated  at  tirr»~  lew!  n 


Equations  (6)  can  be  birt|.**r  «i»rpbfi«d  to 
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where 


U’+l  =  U*  +  AU* 

(7c) 
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Equation!  (7)  describe  an  algorithm  that  can  be  used  to  solve  equation  (1)  To  avoid 
large-amplitude  oscillations,  especially  in  shock  wave  solutions,  a  fourth-order  dissipation 
term  is  added  to  the  right  hand  side  of  equation  (7a).  A  fourth-order  dissipation  term  has 
been  chosen  so  that  the  formal  order  of  the  scheme  is  not  disrupted.  The  dissipation  term 
is  of  the  form 

-(w/»)a*V0'u/0*4)i;  *  ~(u,/8)(u,\,  -  4uj«,  +  «u;  -  4U;_,  +  u;.,)  (u) 

The  above  scheme  (equations  (7)  and  (11))  is  stable  for  values  of  u>  in  the  range  0  <  <  1 

•  tending  to  a  linearised  von  Neumann  stability  analysis  (fleam  ar.d  Warming  1976).  This 
implicit  factored  scheme  is  easy  to  implement  when  compared  to  an  algorithm  with  centered 
spatial  differences,  which  usually  involve  the  inversion  of  tridiagonal  and  pentadiagonal 
matrices.  This  brief  summary  highlights  the  simplifications  that  can  be  made  for  nonlinear 
•vslerns  whose  flux  vectors  are  homogeneous  functions  (of  degree  one).  The  most  important 
of  these  simplifications  is  the  splitting  of  the  flux  vectors  into  subvectors  which  correspond 
to  their  characteristic  speeds 

2.2  Two  Dimensions 

Fleam  and  Warming  1976  als<*  applied  this  implicit  factored  scheme  to  the  two- 
dimensional  Eulenan  invisiad  gKsdynamic  equations  The  value  of  this  approa-h  may 
be  seen  from  the  large  number  of  problems  in  the  literature  to  which  it  has  been  applied. 
These  include  transonic  aerodynamics,  i.e  lifting  and  nonliving  of  airfoils  oscillating  in  a 
free  stream  (Hearn  and  Warning  1970  and  Steger  and  Warming  1981),  Couctle  flow,  shock 
boundary-layer  interactions  and  Warming  1978)  and  in  one-dimensional  shock  tube 

flow  (Sieger  and  Warming  1981),  see  Section  3 
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Beam  and  Warming  1978  u«ed  the  Couette  flow  problem  (unateady  flow  between  two 
infinite  adiabatic  parallel  walla)  aa  a  test  prob'em  for  their  two-dimenaional  algorithm.  The 
spatial  accuracy  and  stability  of  the  numerical  algorithm  aa  well  aa  its  boundary  conditions 
were  chosen  aa  the  initial  teat.  The  numerical  solution  when  compared  to  the  analytical 
solution  exhibited  a  good  degree  of  accuracy.  No  numerical  dissipation  was  added  to  the 
numerical  algorithm  for  this  problem.  The  shock  boundary-layer  problem  was  used  as  a 
more  severe  test  for  the  two-dimensional  algorithm.  No  analytical  solution  for  this  problem 
was  available  for  comparison  but  when  compared  to  the  numerical  solution  obtained  by 
other  methods  this  method  was  in  good  agreement. 

Steger  and  Warming  1981  developed  new  explicit  and  implicit  dissipative  finite- 
differences  for  the  one-dimensional  and  two-dimensional  forms  of  the  inviscid  gasdynamic 
equations.  These  different  methods  were  used  to  solve  the  one-dimensional  shock  tube 
problem  (see  Section  3)  and  the  two-dimensional  problem  of  transonic  airfoils  Centered 
spatial  differences  were  also  used  and  compared  to  upwind  differences  (used  for  the  split 
flux  formulation)  The  implicit  upwind  scheme  (the  method  used  in  this  report)  solved  the 
one-dimensional  shock  tube  problem  better  than  the  other  individual  methods.  Although 
the  combined  algorithm  of  an  explicit  upwind  scheme  and  explicit  MacCormack  scheme 
(Steger  and  Warming  1981)  was  an  improvement  on  the  implicit  upwind  scheme.  For  the 
two-dimensional  problem,  the  conventional  implicit  algorithm  using  cen'ered  differences 
was  compared  to  the  implicit  upwind  algorithm.  The  results  obtained  for  the  airfoil  prob¬ 
lem  using  these  two  different  algorithms  were  in  good  agreement  although  there  were  smalt 
oscillations  in  the  implicit  upwind  solution.  These  oscilla  ">ns  were  due  to  the  conserva¬ 
tive  flux  vectors  having  discontinuous  derivatives  i.e.  the  eigenvalues  changing  sign.  By 
adding  blending  term  to  the  eigenvalues  (Steger  and  Warming  1981)  the  escalations  were 
smoothed  out.  Beam  and  Warming  1976  used  ?.  hybrid  scheme  with  a  fourth  order  dissi¬ 
pative  term  (similar  to  (11))  to  solve  this  airfoil  problem.  This  solution  agreed  with  that 
of  Steg'r  and  Warming  1981  and  had  no  oscillations. 
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3.  A  SIMPLE  DJfE-L)LMEN5I0NAL  PROBLEA* 


A  ftringrnt  test  for  non  linear  syitemi  described  by  the  equations  (1)  and  (2)  is  the  one- 
dimensionai  shock  tube  or  Ricmann  problem.  This  is  a  simple  problem  where  a  diaphragm 
separate*  two  regions  of  different  pressures  and  densities.  The  fluid  on  each  side  of  the 
diaphragm  js  an  ideal  gas  which  is  at  rest  at  time  t  ~  0  The  initiri  conditions  f-;r  this 
problem  are  shown  in  Fig  1(a).  The  diaphragm  is  placed  at  z  ~  0  and  p\  >  p?  The 
pressures  are  assumed  to  be  such  that  pt  >  pg. 

Equations  (1)  and  (2)  can  be  solved  exactly  for  this  problem  (T;.kewaki  and  Yabr  1987) 
enabling  a  precise  assessment  of  the  numerical  method.  This  solution  is  presented  in  detail 
in  Appendix  B  and  illustrated  in  Figure  1(b)  for  density,  ](c)  for  pressure,  1(d)  for  total 
energy  and  He)  for  particle  velocity.  Because  of  the  step  changes  in  physical  properties, 
as  shown  in  Figs  1,  the  Riemann  problem  is  considered  a  severe  lest  of  the  capability  of 
numerical  methods  to  solve  the  Euler  equations.  When  the  diaphragm  is  ruptured  (f  >  0, 
see  Fig  1(b)),  an  expansion  propagates  into  the  high- pressure  gas  and  a  shock  wave 
followed  by  a  contact  discontinuity,  piopagates  into  the  low-pressure  gas  Region  2  in  rig 
1(b)  represents  the  area  of  adiabatic  expansion,  region  3  the  contact  discontinuity  and 
regton  4  the  shock  wave.  A  paricuiar  numerical  method's  accuracy  and  performance  can  he 
determined  when  the  exact  solution  and  the  numerical  solution  are  compared,  especially 
at  the  contact  discontinuity  and  the  shock  front. 

The  B  -am  and  sVarming  met„od  waa  (rated  against  this  Rrrmann  problem  tsing  -qua 
tions  (T)  and  three-point  backward  tim. -differencing  (3  =  1  and  (  J).  The  three  exam¬ 

ples  of  implicit  tune-differencing  given  in  Section  2,  the  tra-'roidal  formula,  the  backward 
Eul'r  differencing  and  the  three-point  backward  differencing,  were  all  used  to  solve  the 
Hiemann  problem.  Tnere  were  no  significant  differences  trtwren  the  solutions  given  by 
these  different  time-diffrrenring  schemes,  except  that  the  trapezoidal  formula  was  unstable 
for  t  >  b  HOI.  Hence  only  the  three  point  backward  timc-differencing  was  used  f.,r  this 
section  The  «ch<  nr  in  uncoifdilibnally  stable  for  6  -  1  and  (  ~  \  (Sieger  ami  vV.irmmg 
IT®!) 

The  dissipation  term  given  in  equation  (11),  was  added  to  the  right  hand  side  of 
ticn  (7<i),  with  u>  —  0  5.  One  hundred  points  were  laker,  eadr  side  of  the  diaphragm  with 
din  -  f' 025,  so  that,  r  was  in  the  rang-  2  !>  ■  z  •  2.5  The  mt.o  of  specific  li-ats.  T,  was 
1  -1  and  the  initial  conditions  wore  p  1  /■':  -  )n  li,  v;  ;w  -  ri  ,md  uj  ■  -  us  -  t) 

There  vas  a  fixed  boundary  at  r  -  -t.U  .n-.l  at  t  2V  1  c.  U"-j  -  l.*.  a„d 
bi'n-r  If, vi  w-h-re  <r  refers  to  tlo-  lime  step,  0  a  he  initial  tru-  step,  .72  the  bc-undiirv 
e  2  5  and  .VI  the  b<»unctary  r  2  5 


in 


4.  DISCUSSION  OF  RESULTS 


Using  the  implicit  factored  scheme  and  the  parameters  given  in  Section  3,  the  Riemann 
problem  was  solved  for  different  time  increments  (At  —  0.01,  At  —  0.005,  At  =  0.0G1  ard 
At  —  0.0001)  but  for  the  same  total  lime  (t  —  '.0).  Kach  of  the  solutions  Fig.  2  5  ca.« 
be  compared  to  the  exact  solution  shown  in  Fig.  Ifb)-l(e).  The  parameters  and  solution 
values  for  the  exact  solution  are  given  in  Table  i. 

The  density  profile  after  100  time  step*  l At  -  0  01)  is  shown  in  Fig.  2(a).  This 
profile  follows  the  general  form  of  the  exact  solution  (Fig  1(b))  except  for  a  non-physical 
oscillation  at  the  shock  front  z  =  z\  and  the  smearing  of  the  contact  discontinuity  and  the 
shock  front.  Although  both  the  contact  surface  and  the  shock  front  have  both  been  smeared, 
this  numerical  scheme  has  resolved  the  shock  front  better  i.e.  there  is  less  dispersion  at 
z  ~  £4  than  at  z  =  zj  By  adding  more  dissipation  (i.e.  increasing  the  value  of  u;  in 
equation  (11))  the  non-physical  oscillation  or  the  Gibbs  type  error  effect  can  be  reduced. 
To  reduce  the  amplitude  of  the  oscillation  more  dissipation  was  added  in  ♦his  way,  but  this 
was  unsuccessful.  The  dissipation  term  is  two  orders  of  magnitude  smaller  than  the  overall 
scheme,  hence  by  increasing  u>  there  should  be  no  significant  chi  in  the  solution  as  was 
observed.  Recall  from  Section  2,  the  addition  of  a  dissipation  *  m  should  not  effect  'he 
overall  stability  of  the  scheme  and  therefore  hta  to  be  ‘his  small.  Hence  to  improve  the 
overall  accuracy  of  the  solution,  which  would  also  rrdn  :e  the  amplitude  of  the  oscillation, 
the  time  step  needs  to  be  made  smaller  i.e.  increas  ;..e  number  of  lime  steps  Before 
considering  a  smaller  time  step  it  i»  worthwhile  to  look  at  this  non- physical  oscillation 
in  the  corresponding  profiles  of  pressure,  energy  and  velocity  in  Fig.  2(b),  fig  2(c)  and 
Fig.  2(d),  respectively  Note  that  the  oscillations  cedi'-  at  the  sharp  change  in  pressure,  but 
not  when  the  density  changes  with  constant  pressure.  The  oscillation  is  more  pronounced 
in  the  pressure  and  energy  profiles  and  is  at  it?  extreme  in  the  velocity  profile  where  it  is 
producing  negative  velocities. 

The  time  step  has  been  reduced  by  a  factor  of  a  half  (A*  =  0  i/J5)  in  Fig.  3.  The 
amplitude  cf  the  oscillation  haa  been  significantly  reduced  with  only  a  small  undershoot, 
consisting  of  three  grid  points,  remaining.  Also  note  that  the  boundaries  z  =  and  z  -  z 4 
have  not  suffered  any  further  dispersion.  Inducing  the  time  s*ep  again  (At  ~  0.001),  see 
Fig.  4,  reduces  the  amplitude  of  the  remaining  undershoot  at  the  shock  fro'»t  but  the 
solution  ha*  become  more  diffuse  with  the  contact  surface  and  the  shock  front  further 
dispersed.  This  has  caused  the  solution  to  become  less  nccurate  when  compared  to  the 
exact  solution.  Thi 6  situation  worsens  when  the  time  step  is  reduced  again  (At  -  0,0001), 
see  Fig  5,  where  both  the  contact  surface  and  the  shock  fr<»n*  have  b*.c"  severe’v  dispersed 
No  oscillations  occur,  but  compared  to  the  exact  solution  in  Fi~.  !»Sj  the  solution  has 
become  inaccurate.  I  h.;  best  balance  brft*ern  the  problems  of  non-ph)  sical  escalations 
and  tevere  diffusion  appears  to  be  in  Fig  J  where  At  0  003  C't^tot.tly  v.nh  different 
initial  conditions  and  parameters  a  lime  step  other  than  0.005  *-  vedd  provide  the  required 
balance  between  diffusion  and  oscillations.  The  best  value  for  At  r  (■  .  objective  compromise 
based  on  diffusion  verses  oscillations  and  c«  only  be  ueiermm'  •  by  .  «a)  and  rn 
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The  Courant*Fredrichs-Lewy  condition  (CFL)  wm  used  to  ‘  stability  of  the 

scheme.  The  CFL  condition  for  this  particular  hyperbolic  scheme  i* 

|A*(U)|^<I,  k,.  1,2,3.  (12) 


where  A*  are  the  eigenvalues  of  the  Jacobian  matrix  .4(U).  Expressing  this  condition 
geometrically:  the  numerical  domain  of  dependence  of  the  scheme  (7)  (with  equation  (11) 
added)  must  contain  the  domain  of  dependence  of  the  differential  equation  (Peyret  and 
Taylor  10S3)  The  solutions  discussed  above  easily  satisfied  this  condition.  This  is  because 
the  numerical  scheme  is  unconditionally  Table  The  CVutanl  number  did  become  close  to 
unity  at  r  «  r4  (the  shock  front),  however.  This  was  the  region  where  the  oscillation* 
occurrd,  hence  venting  that  the  oscillations  were  non*physical  i  e.  they  were  not  part  of 
the  solution  but  an  error  of  the  scheme  employed. 

Ideally  an  infinite  shock  tube  should  he  used  for  the  onc-dunension/tl  shock  tube  prob¬ 
lem  and  therefore  any  end  boundary  effects  could  be  eliminated.  This  is  not  possible  in 
practice  since  the  numerical  scheme  that  has  been  used  here  is  an  alternating  direction 
(AIM)  method  (i.r  equation  (Ta)  is  evaluated  in  the  positive  r  direction  and  then  equation 
(7h)  is  evaluated  in  the  negative  t  direction)  therefore  the  boundaries  at  i  —  —2.5  and 
r  -  2  5  are  fixed  As  can  be  seen  from  the  results  these  boundaries  are  f*r  enough  away 
fri-m  \he  d-.aphrsgni,  hence  there  is  no  boundary  effect  on  the  solution.  With  boundaries 
closer  to  r he  diaphragm  the  stability  of  the  scheme  and  the  numerical  solution  may  have 
been  mlvrt»e!y  ntlected 

*1  he  implicit  factored  scheme  that  lias  I  mi  employ'd  here  crim  'd  be  ure  J  "  her.  the  flux 
vector  F  is  »o;i-lion»ogrnorus.  Itecall  that,  if  K  satisfies  the  homogeneous  properly  Irefcr 
f.,  equation  (1))  and  A  has  n  complete  set  of  hr.carlv  independent  e»g-iiv<-ctors  then  the 
flux  \rctot  F  can  be  *plit  into  two  subvcctai*.  one  subvretor  corresponding  to  the  positive 
eigetuftIi, rs  of  A  and  the  other  the  negative  eigenvalues  These  «ul> vectors  can  then  be 
d nfe i meed  individually  with  the  appropriate  one-smrd  scheme.  If  the  flux  vector  is  non- 
h  un  .geneous  then  it  cannot  be  split  into  suuvectors  and  therefore  an  alternative  algorithm 
iiMUg  cent*  o  !  spatial  ditfererers  has  to  be  used.  As  noted  in  Section  I  Ugoni  funs  using 
upwind  differencing  nave  superior  dissipative  and  dispersive  pn'prrtir;  os  we))  as  being 
more  robust  and  efficient  than  algorithms  umtir  centered  differences.  Therefore,  to  be 
able  k,  solve  in  arbitrary  nonlinear  hvperbdu  >v<tcm  (ti«  Conserv.uion-law  form)  v’lth 
tfiiN  implicit  fav  lop'd  scheme  the  equation  of  slate,  ’.lie  u  uluirar  flux  vector?  must  be 
b-uj*  «u j  J  iiirti-  ns  of  degrrr  one 


K  CONCLUSIONS 


Using  the  one-dimensional  shock  lube  or  Riemann  problem,  the  implicit  factored 
scheme  developed  by  Beam  and  Warming  1976  has  been  tested  for  its  capability  to  solve 
the  Euler  equations.  The  Riemann  problem  proved  to  be  a  severe  test  for  this  scheme.  The 
scheme  performed  well  but  it  had  a  couple  of  drawbacks.  Non-physica)  oscillations  (Gibbs 
error)  and  dispersion  of  the  contact  discontinuity  and  the  shock  front  adversely  affected  the 
solution.  By  removing  the  oscillations  the  solution  became  more  dispersed  and  by  minimiz¬ 
ing  the  dispersion  the  amplitude  of  the  oscillations  increased.  Therefore  the  best  solution 
for  a  problem  (solved  with  this  method)  is  a  subjective  compromise  based  on  diffusion 
versus  oscillations  and  can  only  be  determined  by  trial  and  error.  Besides  these  difficulties 
this  implicit  factored  scheme  was  stable  and  efficient  for  a  wide  range  of  parameters. 

The  main  motivation  for  this  work  was  to  examine  a  paricular  numerical  method  to 
see  if  it  could  successfully  model  shock  waves  and  be  adapted  for  studying  shocks  in  the 
solid  state.  After  finding  a  suitable  numerical  method,  the  simple  one-dimensional  shock 
tube  (Riemann)  problem  presented  in  Section  3  would  be  extended  to  model  shock  waves 
in  solids,  eg.  the  shock  waves  generated  when  a  flyer-plate  impacts  a  solid  material.  This 
requires  different  equations  of  state  to  the  ideal  gas  equation  used  here.  The  ideal  gas  equa¬ 
tion  is  a  special  case  of  the  functional  form  of  the  equation  of  state  (given  in  equation  (2)), 
which  this  method  is  suitable  for.  If  the  equations  of  stale  for  solids  can  be  constructed  in 
this  functional  form  or  if  this  functional  form  is  used  over  restricted  *.*ngcs  of  piessure  and 
density  then  the  implicit  factored  scheme  can  be  used.  If  the  equations  of  state  for  solids  are 
not  in  this  form  then  an  alternative  numerical  scheme  must  be  found.  Elastic-plastic  con¬ 
ditions  must  be  incorporated  into  the  problem  and  usually  the  addition  of  suitable  source 
or  sink  terms  to  equation  (1).  Extensions  to  two  dimensions  and  the  possible  incorporation 
of  initiation  and  detonation  of  explosives  would  follow.  These  problems  would  obviously 
require  a  fairly  robust  and  versatile  numerical  scheme.  The  implicit  factored  scheme  used 
here  performs  well  for  the  sample  problem  but  the  restrictions  or.  the  boundary  conditions 
and  the  choice*  of  equations  of  state  arc  too  restricive  for  our  purposes. 

Though  not  yet  examined  »r  detail  by  us  other  methods,  such  as  flux-corrected  transport 
(FCT)  (Boris  and  Book  1973,  Book,  Boris  And  Ham  1975  and  Boris  and  Book  1976a)  «nd 
Zalesak  1979  multidimensional  FCT  appear  to  provide  the  accuracy  and  versatility  ne  *ded 
for  these  problems.  Results  published  elsewhere  (Boris  and  Book  1973,  Book,  B*'r»s  and 
Bain  1975,  Boris  and  Book  1976a,  1976b,  Zalrsak  1079,  1 98 1  and  Sod  197b)  suggest  that 
they  arc  highly  accurate  and  reasonably  easy  to  use  An  excellent  example  of  this  it  in 
Zalrsak  1981.  The  one-dimensional  shock  tul  *  (Riemann)  problem  is  solved  using  a  J  t'T 
method  with  similar  parameters  to  the  problem  solved  here.  The  EC  T  method  solves  the 
shock  tube  problem  extremely  writ  The  solution  is  highly  accurate  and  has  superior  result* 
when  compared  to  the  method  turd  here  and  other  methods  elsewhere  (Sod  1978)  Further 
investigations  arc  taking  place  ah  ng  these  lines. 
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TABLE  1 


Exact  solution  values  for  the  Riemanrt  problem  (  see  Section  3  ).  The  shock  velocity 
TJ  is  1.90,  see  Appendix  B-  Note  that  across  the  contact  discontinuity  x  —  rj,  pj  =•  p<  and 
u$  -=  u«.  The  initial  conditions  are  p\  =  10.0,  ps  =  1-0,  pi  =  10.0,  p$  =  1.0,  while  7  —  1.4 
throughout. 


Region  i 

Density 

pi 

Pressure 

Pi 

Total 

Energy 

E. 

Particle 

Velocity 

u* 

— 

Boundary 

1 

10.00 

10.00 

25  00 

0.00 

-1.32 

2 

. 

. 

. 

* 

-0.02 

3 

4.08 

2.85 

9  05 

0.97 

4 

2.04 

2  85 

8.09 

097 

1.90 

5 

1  00 

1.00 

2.50 

0  00 

' 
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Riemann  Problem 


Figure  1(b)  Density  t  rofile  after  thr  diaphragm  Iim  beer.  ruptured  Also  u  the  foim  of 
the  exact  solution,  see  appendix  B. 


(c>  t  >  0 

Figure  1(c)  Same  aj  Fig.  1(b),  except  preuure  profile. 


(d)  t  >  0 


Figure  1(d)  Same  u  Fig.  1(b),  except  total  energy  per  unit  volume  profile 


Figure  1(e)  Same  Fig.  He),  except  particle  velocity  profile 


Figure  2(a)  Riemann  problem  solution  using  implicit  factored  sch< 
i  ~  1  0  with  At  =  0.01,  Ax  =  0.025  and  100  points  each  side  of  the 


figure  2(d)  Same  u  Fig.  2(a),  except  particle  velocity  profile. 


Figure  4  Same  u  Fig.  2(a),  except  At  ~  0.001. 


In  this  appendix,  the  concise  form  of  the  implicit  factored  s -rheme  presented  in  Se  :tion  2 
is  expanded. 

The  one-dimensional  form  of  Euler's  equations,  namely,  the  continuity  equation,  con¬ 
servation  *:(  mass  and  conservation  of  energy  are  represented  m  equations  Ml;,  (A'J)  and 
(/t3),  respectively  given  by 

^  +  ~(pu)  =  0  (41) 

'(£-£)+2- 

'U  u£)  +  l:fpu)  =  0  M3) 


where  e  is  the  internal  energy  per  unit  mass,  p  -  (7  —  l)pe  from  the  ideal  gas  equation  of 
state  ana  7  the  ratio  of  specific  heats. 


These  equations  can  be  written  in  vector  form  as 

<21! ,  2E1H2  -  0 

9t  dr 


(-11/ 


where 


and  F.  ii  the  total  energy  per  unit  volume  and  p  -  (7  -  1)(E  -  \pn7) 

Tl.r  K\*!:r;rn  equations  of  gasdyn&mics  have  the  property  that  the  nonlinear  fuii'-tion 
F(U)  is  a  homogenous  function  <■{  degree  one  .n  U,  sec  Beam  and  Warming  l£7b 

Therefore, 

F(U)  ----  A(V)V  (  \ r> I 


w  here 


A(V)  = 


0  F 

d'v 


A“  - 


i»  the  Jacobian  matrix.  Hence  equation  (A4)  can  be  written  a> 


(A7) 


At  noted  in  Section  2  there  ate  numerout  example*  of  implicit  time  differencing  A 
•mgle-atep  temporal  «cheme  for  advancing  the  tolution  of  equation  (.At)  (Beam  and  Warm¬ 
ing  1978)  it 


U" 


dAt  t/fHV  <0' V 

ntKft)  +(«) 


dAt 

f  - 


At 


1  +  (  '  Ot 


(?H)* 

V  fit  > 


1  +{ 


'] 

(U--U*-‘)d  0((f-  J-OAl*  +  Ats).  (,1ft) 


Here  U(l)  =  U(nAf)  =  U*  and  At  it  the  ditctete  time  increment,  d  and  (  define  the 
particular  time-difference  approximation  (tee  Section  2).  Substituting  equation  (,U)  into 
(.-18)  givrt 


U'tl  ir 


Udr\' 

,f)F. M  *1 

At  ,0Fv' 

[(fix)  + 

(fix )  1“ 

J-(U"  -  U-')  +  tV(«-  j  -OAt5  4 


At3]. 


Thu  equation  it  nonlinear  due  to  the  term  involving  F" ' 1  -  F(U"  +  1 )  In  order  to  solve 
for  tr"  in  equation  (T9)  the  nonlinearity  mint  lie  removed.  A  local  Taylor  expand!  ,n 
about  U*  give# 

r‘,!sr  +  Q  V'H  O(M')  (,tl0) 

Mil'*tilutin&  in  equation  (.16) 

F”“  -  F-  t  .■r(V'n  -  U")  ,  <>(A<5)  (All) 


Hrncv. 

r“  .-rv!*’1  .  (),Ar 

I  heref.vre,  equati-m  ,  ,19)  In-coirrc 

U’1*'  V 


d.\t  I  ft 


’  k()t  j  1  ♦  V  > 


(U2} 


,V,<ir  {  -O-V1  .  A<’j  c  u:u 


a  -  : 


Rearranging  equation  (-413)  girt* 
Mt  6  ' 


ft .  t^L.  /im] _ +  X-au”'1 

lI  +  1  +  (d*  ~  i  +  <W  +i+{AU 


(-414) 


where 


AU*  =  U*+1  -  U\ 


(-415) 


Thi»  equation  is  a  noniterative,  setond-ordet  time-accurate  solution  for  equation  (/l-t). 

The  Jacobian  matrix  -4  ran  be  simply  calculated  to  be 

/  0  1  ON 


4  = 


b=j*' 


(3  -  7)u  7-1 

l(T-l)»’  -  ^  ■^-^(T-l )«’  7»  ) 


(-410) 


The  eigenvalue*  of  A  are 


At  =  a,  Aj  =  u  +  c,  Aj  =  u  -  c 


(417) 


where  e  =  (7^)^  is  the  local  speed  of  sound. 

Since  F  is  a  homogeneous  function  (see  equation  (45)  and  Steger  and  Warming  1981), 
F  van  be  split  into  two  paits 

F  =  F+  4  F~  (418) 

where  F  7  corresponds  to  the  positive  eigenvalues  of  4  and  F~  the  negative  eigenvalues  of 
4 

Any  eigenvalue  A|  can  be  expressed  as 

Aj  =  V+Af  (419) 


where 


,+  A,  +  |A(| 

A  =  g  . 


a- 


(420) 


so  if  A|  >  0,  then  A,+  >  0,  A,”  -  0  and  il  A;  v.  0  then  A,'  0,  =  0. 

For  the  system  described  by  (48)  to  be  hyperbolic  there  must  exist  a  similarity  trans 
formation  such  that 


/  Ai 


V  0 


Am  / 


A- 3 


Q-'AQ  -  A 


(421) 


where  A  is  a  diago-.ui  matrix  end  the  eigenvalue*  of  A  ere  read. 

Hence  the  diagonal  matrix  A  can  be  expressed  as 

A  =  A+  +  A"  (A??) 

where  A+  and  A"  have  diagonal  elementa  A(+and  A,~,  reapectively 
Therefore, 

A**QA+Q-'  A-  =  Q\-Q~' 

Ff  =  A*U  f-  =  A'U 

and 

A  =  A1-  +  A~ . 

If  g-  ia  replaced  by  the  backward  difference  operator 

V,  =  U,  -  U,.,  (/126) 


(.423) 

(•424) 

(.425) 


then  (d8)  ia  atahle  if  and  only  if  the  eigenvaluea  of  A  are  all  positive. 

Similarly,  if  gt  ia  replaced  bv  the  forward  difference  operator 

A,  U,M~U,  ( .427) 

then  (.48)  is  stable  if  and  only  if  the  eigenvalues  of  .4  are  all  negative. 

A  noniterative  implicit  finite  difference  a  rheme  (Beam  and  Warming  1 976),  for  a  one. 
dimensional  system  of  conservation  laws,  with  the  use  of  split  flux  vedota  (.423)  and  one¬ 
sided  spatial  difference  approximations  ia 


■(^„TVp)(<:F.4r  +  (.,F;r) .  jd-jiu;-'  (.ran 


where 


And 


elK.  2F,  -  4Fj- ,  +  F, - j 
■\'K.  -  3F,  v  :p„,-  f,(J 


By  redefining  some  of  these  f.vlor*.  that  is. 

BM 

(>i  Ar(l  -.  () 


A-4 


At 

“*  =  2Ar(l  +  {) 


•mi 


as  = 


( 

l  +  < 


the  *cheme  (428)  can  b«  implemented  ai 

(i  +  «,v,/t;r)Au;  =  -a,(«tF;r  +  +  ojau;-1 


(429o) 


(i+a^.^rjAUJraU;  (/1 296) 

o';  r>  =  U*  4  Au;.  (,429c) 

Further  «;r..plifing  of  (429)  gives 

i/1-Au;  =  o,.4;.,i"au;.,  -o,(^f;i*  +  s'Fjr)  +  «sau;-'  <4.io«) 

i’au;  -  au;  -  o,.i;ui'Au;M  (43i>6) 

U"+1  =  u;  +  AU;.  (,430c) 

where 

i,;r  -  r  +  o,,i;r  . .  (4.10 

nn<i 

L;p  -  I  -  a,.\;p  (.432) 


The  scheme  elrscftbed  by  (.-130)  -  (.132)  was  useti  \o  5<*lve  ihe  one -dimensional  problem 
in  Sertior  3. 


A -5 


fl.  APPENDIX  B 


In  this  appendix,  the  exact  lolution  for  the  problem  described  in  Section  3,  al»o  known 
u  the  Riemann  problem,  ii  derived.  The  derivation  follow*  that  of  Tnkewaki  and  Yahe 
19d7,  except  a  relation  wa*  obtained  to  find  the  (hock  velocity  V  instead  of  the  Marh 
number  Af.  The  initial  condition*,  namely  a  diaphragm  placed  at  r  =  0  (between  the 
region*  1  and  5),  are  shown  in  Fig  1(a)  This  profile  i*  modified  into  that  in  Fig  1(b),  when 
t  >  0. 

The  Rankine-Hugoniot  conditions  (Haye*  1973)  aero**  the  (hock  front  x  =  r4  are 


Pi  . 

Pi  v 


(in) 


and 


n-  n  -  PiUin 


Ei  ~  F-i  =  -(p*  +  r*)(»>!  -  e4) 


(B1) 

(in) 


where 


For  a  perfect  ga» 


therefore  expiation  (bi)  become* 


p*v*  -  p5l'5 


7  ~ 
2 


(Ps  +  PaXu*  -  vt). 


(n  4) 


From  equation  (#1), 

«*  --1/(1  - Pi) 

pi 


(in) 


Substituting  equation  (775)  into  (lrl)  gives 

pi  -  rs  ~  b'7(pi  -  pi) 

Pi 


(/7fl) 


From  equation  (Bi), 


P>  (7  *  Dpi  -  (7  -  1  )pr 

Pi  (7  *  Dps  -  (7  -  Dp* 


(B7) 


B-i 


By  aubatltuting  thii  equation  into  (Bfl),  an  expretaion  for  p\  in  termi  of  the  known  quan¬ 
tity  in  region  5  ii  obtained, 


_ e\V'(t  +  i) _ 

P<  *  27Pi  +P»(,!( ~~  O' 


(//«) 


Substituting  equation  (D8)  into  (B 7)  girea, 

ip*U*  -  (7  -  l)p« 

p,“ — (7+Tj — 

und  Rubfttituting  equation  into  (fl$)  gives, 

/rpi) 

*  Pat/(7  +  0 


(DO) 


(D 10) 


The  preaaure  and  the  velocity  ahould  be  continu  >u»  at  the  contact  aurface  * 
Hence, 


U|  «  u« 


(All) 


and 

Pi^Pi  (H12) 

The  particle  velocity  u  i<  ataurned  to  lie  a  aiiiRle- valuer!  function  oflbe  denaily  p  in  the 
adiabatic  expansion  region  (region  2)  Therefore  eqnationa  (!)  and  (2)  reduce  to 


Op 

Ot 


+ 


7P\fj  Op 
p  ’  \  0 


(D 13) 


if  the  function  u(p)  aaliafiea  a  relation 


du  _  _  I  ( 7P\  i 

Jp  p\  p  ) 


(DU) 


From  equation  (W13),  the  veloritiea  of  the  rarefaction  front  1  *  ■  i  |  )  and  the  boundary 

at  r  x  X}  are  -(7Pi/pi)^  and  uj  -  (jpj/p j)^,  reapectively  t ntc-p.rtat the  Inal  equation 
yivea 


(D 15) 


Since  the  adiabatic  rcatiun 


Pi 


r  i\t 


hold*,  ut  i*  related  to  pi  u 


(B 17) 


Finally  by  *ub*tit-iting  equation*  (/? 6),  (OlO),  (All)  and  (012)  into  (B 17),  give*  the 
relation  which  determine*  the  ihock  velocity 


(y  -  1  )(/>fU3  -  ipi)f  pi  \>  m  _  (IpiU*  -  (?  -  I )p* N  V 

P|l/(7 -e  I)  V-r|»i  /  "  V  (7  +  l)pi  / 


(018) 


When  equation  (018)  tea*  *o)ved  u*lng  (he  tame  parameter*  in  Fig  5  (Table  1)  of 
Takewakl  and  Yabe  1987,  the  tliock  veiocitte*  agreed.  The  .lumerical  value*  of  the  exact 
•nlution  for  the  problem  in  Section  3  are  given  in  Table  1. 
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